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ABSTRACT 


A  seismic  refracCion  survey  was  conducted  at  Reydarf jordur ,  Iceland  as 


part  of  the  Iceland  Research  Drilling  Project  in  July,  1978.  An  array  of  15 
seismometers  was  used  to  measure  apparent  velocities  from  25  explosions  along 
a  ,23  km  refraction  line.  The  experiment  yielded  data  used  to  determine  the 
P  and  S-wave  velocity  structures  beneath  the  drill  site  and  to  obtain  bounds 
on  the  depth  to  seismic  "layer  3".  The  observed  apparent  velocities  vary 
continuously  as  a  function  of  range  with  a  marked  increase  in  velocity  at  a 
range  of  16  km.  This  corresponds  to  an  observed  cross-over  in  the  travel  time 
data.  The  measurement  of  apparent  velocities  has  yielded  data  of  the  form 
X(p),  the  distance  to  a  measured  ray  parameter,  and  T(p),  the  delay  time. 

Joint  inversion  of  the  X(p)  and  T(p)  data  has  been  performed  by  linear 
programming.  Layers  of  constant  slowness  gradient  are  used  as  the  basis  for 
the  inversion  process.  Estimates  of  error  in  the  data  are  used  to  produce 
extremal  bounds  on  the  family  of  earth  models  that  satisfy  the  data. 

The  resulting  envelope  of  P-wave  velocity  models  indicates  a  continuous 
variation  of  velocity  with  depth.  Two  regions  of  high  velocity  gradient  are 


observed,  one  near  the  surface  and  the  other  at  depths  between  2.6  and  3.9  km 
% 

associated  with  "layer  3"  in  Iceland.  These  results  suggest  that  modeling  of 
Iceland's  upper  crust  with  a  few  thick  homogeneous  constant  velocity  layers  is 


an  inadequate  representation  of  the  velocity  structure. 
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INTRODUCTION 

The  primary  o  j  ■  Live  of  the  seismic  experiment  at  Reydarf jordur  was  to 
determine  the  P  an  S-wave  velocity  structure  underlaying  the  Iceland  Research 
Drilling  Project  s  te.  The  refraction  survey  was  laid  out  to  allow  determina¬ 
tion  of  a  velocity  model  for  the  stratigraphy  sampled  by  the  1.9  km  deep 
borehole  and  to  allow  comparison  of  the  results  of  field  seismic  methods  with 
those  of  lab  analysis  of  rock  physical  properties.  In  situ  measurements  of 
seismic  velocities  provide  us  with  a  more  integrated  velocity  vs.  depth  model, 
averaged  over  local  heterogeneities  as  opposed  to  the  discrete  sampling  done 
in  laboratory  analyses  of  rock  velocities. 

The  refraction  experiment  was  designed  to  produce  a  bounded  estimate  of  the 
depth  to  "layer  3"  velocities  in  the  vicinity  of  the  drill  site.  Palmason 
(1963,  1970)  suggested  that  "layer  3"  velocities  may  be  within  3  km  of  the 
surface  in  northeastern  Iceland  and  possibly  shallower  beneath  tertiary 
volcanic  centers.  The  close  proximity  of  Thingmuli,  a  Tertiary  volcano 
approximately  5  km  west  of  the  drill  site,  suggests  that  this  seismic  layer 
may  be  close  to  the  surface  along  the  refraction  profile.  This  possibility 
was  investigated  by  generating  and  recording  waves  refracted  by  this  layer. 

A  line  with  shot  to  receiver  separation  up  to  23  km  was  used  to  sample  "layer 
3"  velocities.  Array  recording  of  closely  spaced  shots  was  used  to  provide 
resolution  for  measuring  velocity  gradients  as  a  function  of  depth. 

Finally,  this  paper  investigates  the  appropriateness  of  homogeneous 
layered  Interpretations  of  refraction  data.  The  use  of  thick  planar  layers  of 
constant  velocity  for  the  purpose  of  modeling  travel  time  data  has  been 
prevalent  for  many  years  in  seismic  refraction  interpretation  although  there 
are  many  reasons  to  believe  that  velocity  varies  continually  with  depth.  A 
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thick  layered  crustal  model  generally  represents  an  overconstrained  solution 
to  the  Inverse  problem  and  is  an  artifact  of  the  modeling  technique  and  sparsity 
of  data.  Therefore,  the  Inversion  applied  to  the  data  here  allows  the 
determination  of  a  velocity  function  that  varies  continuously  as  a  function 
of  depth. 

Experiment  Description 

The  refraction  survey  at  Reydarf jordur  utilized  a  densely-spaced  array  of 
seismometers  and  an  unreversed  line  of  shots.  The  use  of  a  fixed  array 
minimizes  error  in  the  measurement  of  apparent  velocities  because  individual 
station  corrections  can  be  easily  determined.  Measurement  of  the  apparent 
velocities  across  an  array  is  equivalent  to  the  direct  measurement  of  p,  the 
ray  parameter.  The  resolving  power  of  the  array  is  dependent  on  the  array 
length,  number  of  sensors  and  frequency  of  seismic  energy  propagating  across 
it.  A  shorter  array  decreases  resolution  in  measured  p  while  a  longer  array 
smooths  out  the  variability  of  the  measured  ray  parameters.  In  refraction 
surveys  where  a  closely  spaced  array  is  not  used,  p  *  dT/dX  (the  slope  of  the 
time-distance  curve)  is  often  approximated  by  fitting  a  smooth  curve  to 
sparse  T,  X  data.  This  interpolation  approach  may  introduce  artifacts  into 
the  interpretation.  Therefore,  we  have  used  the  array  measurements  of  ray 
parameter  for  the  determination  of  velocity  structure. 

The  refraction  survey  was  conducted  on  a  line  west  to  east  with  shots 
ranging  from  1.5  km  to  20  km  east  of  the  drill  site  (Figure  1).  A  15-channel 
digital  recording  array  was  used  in  two  configurations  to  yield  T,  X  data  and 
phase  velocities.  In  the  first  configuration,  the  seismometers  were  placed 
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in  a  closely  spaced  array  centered  about  the  drill  site  for  a  total  spread  of 
2.4  km.  Eleven  shots  were  fired  between  1.5  km  and  8  km  east  of  the  drill  site. 
These  explosions  were  detonated  in  shallow  water  to  enhance  acoustic  coupling. 


Following  these  shots,  nine  of  the  seismometers  were  moved  farther  vest 
of  the  drill  site  to  create  an  array  of  6  km  spread.  These  were  located  on  the 
flank  of  Thingmuli  volcano  at  considerably  higher  elevation  than  the  other 
stations.  This  move  enabled  recording  of  refracted  rays  whose  bottoming 
points  were  close  to  the  drill  site  to  sample  the  velocities  beneath.  In  this 
part  of  the  experiment,  shots  were  located  in  the  fjord  from  7  km  to  20  km  east 
of  the  drilling  project  in  10-20  meters  of  water. 

Joint  Inversion  of  t (p)  and  X(p)  data 

The  analysis  of  the  Reydarf jordur  data  is  based  on  the  measurements  of 
apparent  velocity  with  the  array.  The  data  have  been  reduced  to  estimates  of 
X(p),  the  distance  to  a  given  ray  parameter  and  t(p),  the  delay  time  function. 
Given  an  Infinite  number  of  X(p),  T(x),  or  t  (p)  data  containing  no  observational 
errors,  it  would  be  possible  (in  the  absence  of  low  velocity  zones)  to  find  a 
unique  earth  model  satisfying  the  data  exactly.  Vith  data  sampled  at  discrete 
points  it  is  only  possible  to  place  bounds  on  a  family  of  solutions.  The 
existence  of  observational  errors  in  the  data  further  increases  the  range  of 
possible  models  that  fit  the  data.  Numerous  inversion  methods  have  been 
devised  to  treat  bounds  on  either  X(p)  data  (e.g.,  McMechan  and  Wiggins,  1972) 
or  t (p)  data  (e.g.,  Bessonova  et  al.,  1974,  1976)  and  to  produce  estimates  on 
the  range  of  earth  models  that  satisfy  the  constraints  of  the  data.  However, 
because  bounds  on  the  X(p)  and  r (p)  data  represent  Independent  information, 
it  is  desirable  to  utilize  them  together  as  simultaneous  constraints  in  the 
inversion  process. 

The  inversion  applied  here  utilizes  a  scheme  for  the  joint  treatment  of 
X(p)  and  t(p)  data  as  formulated  by  Orcutt  (1980)  following  the  linear 
programming  approach  of  Garmany  et  al.  (1979).  In  this  approach  a  stack  of 
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layers  of  const  n  radient  of  slowness  is  used  to  approximate  a  continuous 
velocity  vs.  de,  t  -iation.  The  functions  X(p)  and  t (p)  are  represented 
as  linear  functio  .>!  jf  s(p)  (Dorman,  1979)  and  the  seismic  inversion  solved 
using  linear  prof  :.mming.  This  produces  extremal  bounds  on  the  depth  to  a 
given  value  of  p  subject  to  the  constraints  of  the  data  and  the  physical 
^'realizability”  of  the  solution  (Garmany  et  al.,  1979). 

Determination  of  the  X(p)  and  t (p)  data  was  accomplished  by  the 
construction  of  record  section  plots  of  all  seismograms  recorded  for  a  single 
shot  (e.g..  Figure  2).  The  regression  line  T  *  r  +  pX  was  fit  to  the  first 
arrival  picks  (T,  X)  using  the  least  squares  criterion.  The  time  intercept 
of  the  line  is  our  estimate  of  tau  and  the  slope  of  the  line  is  our  estimate 
of  p  over  the  array  interval.  The  value  of  X  was  chosen  to  be  equal  to  the 
mean  shot-to-receiver  distance  to  produce  the  X(p)  estimates.  Some  secondary 
arrivals  were  inferred  from  the  record  sections  to  provide  information  on  the 
extent  of  the  triplication. 

Because  the  instruments  were  used  as  a  closely-spaced  array  for  only  half 
of  the  shots,  an  alternate  method  was  used  for  determination  of  X(p)  for  the 
configuration  where  the  receivers  were  more  widely  spaced.  A  second  order 
polynomial  was  fitted  to  the  (T,  X)  data  using  the  least  squares  criterion. 

The  ray  parameter,  p,  was  computed  for  the  mean  of  X  and  compared  to  p 
computed  by  linear  regression.  The  values  of  p  computed  by  polynomial 
regression  did  not  vary  significantly  from  those  computed  by  linear  regression 
so  the  latter  values  were  used  in  the  modeling  and  inversion  processes. 

To  facilitate  the  placement  of  bounds  on  the  X(p)  data,  confidence  limits 
for  p  were  computed  at  the  95Z  probability  level  assuming  a  normal  distribution 
of  the  (T,  X)  data  about  the  line  T  -  t  +  pX.  These  limits  were  used 
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together  with  th  .*  ■  .it  ter  in  the  X(p)  data  to  infer  upper  and  lower  bounds  for 
the  envelope  of  :1  ;sible  X(p)  values  as  shown  in  Figure  3.  Values  of  x  (p) 
wore  determined  *1  ilarly  where  x  is  the  time  intercept  of  a  regression  line 
fitted  to  the  (T,  )  data.  The  values  of  x (p)  for  the  secondary  arrivals  was 

Inferred  from  rect  rd  section  plots  of  the  seismograms  and  these  points  were 
added  to  the  x(p)  data  to  represent  the  retrograde  branches  of  the  travel 
time  data.  The  scatter  in  all  of  the  x(p)  data  provided  a  final  constraint 
on  the  placement  of  bounds  on  the  x(p)  curve  as  shown  in  Figure  4. 

Compressional  Wave  Velocity  Profile 

The  extremal  slowness-depth  model  for  the  compressional  waves  was  obtained 
by  the  joint  inversion  of  the  X(p)  and  x(p)  data  bounds  contained  in  Figures  3 
and  4.  These  bounds  on  depth  to  a  given  slowness  are  plotted  in  Figure  5.  Two 
zones  of  high  slowness  gradient  appear  in  the  model,  one  in  the  top  .5  km  of 
crust  and  the  other  between  depths  of  2.6  and  3.9  km.  The  former  can  be  attri¬ 
buted  to  surface  weathering  and  porosity  variation.  Because  porosity  decreases 
with  depth  due  to  greater  overburden  pressure,  the  observed  velocities  increase 
as  a  function  of  depth.  The  second  zone  of  high  slowness  gradient  can  be 
identified  as  the  "layer  3"  boundary.  Although  this  boundary  has  been  commonly 
modeled  as  an  abrupt  Increase  in  velocities,  it  appears  that  a  transitional 
region  of  Increase  occurs  over  nearly  1  km  in  depth.  This  contention  is  based 
on  the  small  X-extent  of  triplication  in  the  travel  time  curves  and  the  lack  of 
obvious  reflections  and  secondary  arrivals  from  this  region. 

The  structure  determined  here  does  not  contain  the  features  of  homogeneous 
layers  often  modeled  in  refraction  surveys.  Instead,  slowness  decreases 
continually  as  a  function  of  depth  as  evidenced  by  decreasing  values  of  ray 
parameter  with  increasing  distance.  This  results  from  the' use  of  analysis  that 
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requires  a  continual  velocity  depth  variation  and  not  from  the  technique  of 
modeling  a  finite  number  of  layers  with  discrete  values  of  velocity. 


Shear  Wave  Velocity  Profile 

Because  of  the  sparsity  of  good  S-wave  data,  it  has  been  necessary  to  fit 

a  simple  layered  solution  to  the  S-wave  data.  The  travel  times  of  specific 

phases  of  the  emergent  S-waves  were  measured  and  a  standard  "layer  cake" 

inversion  was  performed.  This  has  produced  a  three-layered  model  of  constant 

velocity  layers  for  the  S-wave  structure  (Figure  6) . 

Knowledge  of  V  and  V  are  aids  in  determining  the  physical  properties  of 
s  p 

rocks  at  depth.  The  values  of  determined  here  agree  fairly  well  with  those 
of  Pdlmason  (1971),  although  the  "layer  3"  velocities  here  are  slightly  higher. 


Because  of  the  uncertainty  associated  with  the  measurement  of  V  and  the  continual 


variation  of  V^,  the  computed  values  of  Poisson's  ratio  exhibit  a  range  of  values, 
from  .26  to  .30.  Palmason's  values  for  layers  1,  2  and  3  are  all  approximately 
.27  with  standard  deviations  of  .02  due  to  a  much  greater  number  of  samples. 


Two-Dimensional  Structure 

The  refraction  experiment  at  Reydarf jihrdur  yielded  Information  that 

% 

can  only  be  used  to  estimate  a  single  dimension  of  velocity  variation.  However, 
one  feature  of  the  data  suggests  the  possibility  of  a  dip  of  "layer  3"  in  this 
area.  The  apparent  velocities  measured  for  waves  refracted  by  this  layer  are 
consistently  higher  than  values  reported  elsewhere  in  Iceland.  The  apparent 
velocities  reported  here  vary  from  6.6  to  7.1  km/sec,  increasing  with  range 
while  Pilmason  (1971)  found  velocities  of  6.35  km/sec  with  a  standard  deviation 
of  .21  km/sec  based  on  72  independent  measurements.  The  higher  velocities 
observed  at  Reydarf jbrdur  could  be  attributed  to  an  eastward  dipping  interface 
beneath  the  west  end  of  the  refraction  line.  Only  a  moderate  degree  of  dip 
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(roughly  3  degrees)  is  necessary  to  explain  the  increase  in  apparent  velocities. 
This  supports  Pdlmason's  view  that  "layer  3"  in  Iceland  may  be  nearer  the 
surface  beneath  Thingmuli  and  other  major  Tertiary  volcanic  centers. 

The  approach  taken  in  this  analysis  has  assumed  only  one  dimension  of 
velocity  variation.  Although  there  is  some  evidence  of  a  more  complex  structure, 
the  lack  of  a  reversed  refraction  profile  precludes  a  two-dimensional  modeling. 
The  presence  of  dipping  or  nonplanar  velocity  structure  may  slightly  alter  the 
depths  to  a  given  velocity  but  the  observed  gradients  will  be  the  same. 

A  Shallow  Low  Velocity  Zone 

A  casual  inspection  of  Figure  7  shows  that  the  travel  times  obtained  from 
the  velocity  depth  model  derived  from  apparent  velocities  across  the  array  do 
not  fit  the  time-distance  data  particularly  well,  and,  in  addition,  there  is  a 
rapid  attenuation  of  the  first  prograde  branch  beyond  10  km  distance.  This 
attenuation  does  not  seem  to  be  correlated  with  the  experimental  procedure 
(movement  of  the  seismometers,  for  instance)  and  it  has  a  continuous  behavior 
so  that  one  would  presume  that  it  is  real.  It  is  conceivable  that  a  form  of 
lateral  heterogeneity  might  be  responsible  for  these  observations,  but  the  gross 
features  can  also  be  modeled  with  a  laterally  homogeneous  structure  with  a 
shallow  and  rather  extreme  velocity  inversion. 

The  reflectivity  method  was  used  to  calculate  synthetic  seismograms  (Fig.  8). 
The  first  model  used  had  a  gradient  zone  that  lay  within  the  extremal  bounds 
down  to  2  km  (the  depth  of  the  borehole)  where  various  low  velocity  zones  were 
inserted.  The  only  model  that  gave  a  reasonably  good  fit  to  the  observed 
first  arrivals  was  one  with  an  abrupt  and  extreme  inversion  (down  to  4  km/s). 
However,  such  a  sharp  change  in  structure  gave  rise  to  prominent  reflections 
at  short  distances.  The  original  data  were  examined  for  sych  reflections,  and 
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it  appears  like.)  at  they  are  present  (Fig.  9).  These  four  pairs 

ve re  used  to  co  is  t  a  Clalraut  envelope  (Bessonva  et  al.t  1974)  from 
which  a  portion  c  t’ ft  delay  time  curve  for  these  reflections  was  inferred. 

The  bounds  to  th:  reflector  were  1.65  to  2.51  km  depth.  It  is  interesting 

to  note  that  a  r  flectivity  run  with  the  low  velocity  zone  at  1.75  km  depth 
yielded  a  much  sharper  attenuation  of  first  arrivals  at  10  km,  and  actually 
looked  more  like  the  data,  but  this  placement  is  not  consistent  with  laboratory 
velocity  measurements  (Christensen  and  Wilkens,  this  volume). 

Although  the  1.9  km  core  does  not  show  any  obvious  lithologic  change  that 
could  be  identified  with  this  reflector  or  velocity  inversion  in  the  lower 
kilometer,  the  other  drilling  results  show  that  there  is  a  sequence  of  highly 
impermeable  layers  beginning  at  about  1.2  km  depth  (Fig.  10).  It  is  conceivable 
that  high  pore  pressure  due  to  water  trapped  under  these  layers  could  lead  to 
an  in  situ  velocity  reduction  which  would  not  be  observed  in  laboratory  velocity- 
measurements. 

N.  Christensen  (personal  communcation)  feels  that  pore  pressure  is  an 

unlikely  means  of  obtaining  the  lowered  velocities,  and  has  suggested  that 

thermal  effects  are  more  significant.  It  is  always  possible  to  argue  that 
% 

refraction  experiments  yield  lateral  averages  of  structures,  and  that  the  lid 
of  the  low  velocity  zone  just  happens  to  be  out  of  reach  of  the  drill  hole. 

Only  a  few  hundred  meters  of  relief  in  the  lid  surface  would  be  required  to  do 
this.  The  particular  velocity  in  the  low  velocity  zone  of  4  km/s  is  conjectural, 
since  it  is  possible  to  trade  off  thickness  versus  velocity,  even,  to  some 
extent,  with  synthetic  seismogram  modeling. 
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DISCUSSION  AND  CONCLUSIONS 

In  the  Initial  design  of  this  experiment  it  was  planned  to  deploy  an  array 
of  seismometers  which  could  be  used  to  obtain  phase  velocities  directly  and 
these  in  turn  used  to  invert  for  velocity  as  a  function  of  depth.  The  rationale 
was  that  this  method  is  insensitive  to  time  perturbations  under  Individual 
shots,  which  would  not  feasibly  be  determined. 

As  the  analysis  of  the  time-distance  data  has  shown  non-uniquness  due  to 
low  velocity  zones  is  a  real  problem  and  a  low  velocity  zone  cannot  be  excluded 
by  the  T-P  data  (Figure  4).  However,  a  comparison  of  the  synthetic  seismograms 
with  the  data  in  Figure  7, in  particular  the  cusd  region  from  12  to  20  km, 
shows  that  the  low  velocity  model  still  is  not  an  ideal  fit  to  the  data. 

If  we  allow  for  the  possibility  of  inhomogeneities  causing  time  delays  under 
the  shot  points,  especially  the  two  shots  at  13  and  15  km,  then  the  low  velocity 
interpretation  might  be  less  attractive. 

With  these  reservations  the  following  can  be  concluded  from  this  study. 

(1)  Modeling  of  the  upper  crust  in  Iceland  with  a  few  homogeneous  constant 

velocity  layers  is  inadequate.  Velocity  varies  continuously  with  depth  apart 

from  the  possibility  of  an  abrupt  decrease  near  2  km  depth.  Velocity  increases 
♦ 

near  the  surface  may  be  attributed  in  part  to  a  decrease  in  porosity  with  depth. 

i 

(2)  There  are  two  regions  of  high  positive  velocity  gradient  in  the  upper 
crust,  one  near  the  surface  and  the  other  at  depth  between  2.6  and  3.9  km. 

The  near  surface  gradient  is  probably  due  to  weathering  and  fracturing.  The 
deeper  region  represents  "layer  3"  and  probably  is  associated  with  a  metaraorphic 
transition. 

(3)  The  contention  that  "layer  3"  in  Iceland  is  locally  shallow  beneath 
Tertiary  volcanic  centers  appears  valid.  However,  the  drill  hole  did  not  reach 
sufficient  depth  to  sample  it. 
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Map  of  refraction  configuration  in  vicinity  of  Reydarf jordur . 

.  Record  section  of  shot  point  8.  All  seismograms  (unfiltered) 
recorded  for  this  shot  at  distance  8.3  km  from  the  array  center. 

.  Envelope  of  P-wave  X(p)  data.  The  solid  lines  represent  smoothed 
bounds  applied  to  the  X(p)  data  and  used  as  constraints  in  the 
inversion  for  velocity  structure. 

Envelope  of  P-wave  r(p)  data.  The  solid  lines  represent  smoothed 
bounds  used  as  constraints  on  the  inversion. 

Bounds  on  depth  found  by  the  joint  inversion  of  X(p)  and  x(p) 
data  (dots).  The  solid  lines  represent  the  inversion  of  the  X(p) 
bounds  alone.  The  dashed  line  represents  the  model  used  for 
determining  the  first  arrival  time  curves  plotted  in  Figure  6. 
Record  section  of  shots  recorded  by  stations  A,  5  and  6.  A 
continuous  velocity  model  was  used  to  determine  the  P-wave 
travel  time  curve.  A  layered  model  was  used  to  determine  the 
S-wave  travel  times. 

Resulting  one-dimensional  velocity  models.  The  S-wave  model  is 
a  siriple  layered  solution.  Palmason's  S  velocities  are  also  shown 
(vertical  bars) .  The  P-wave  model  places  extremal  bounds  on  the 
range  of  possible  depths  to  a  given  velocity. 

Three  sets  of  synthetic  seismograms  derived  from  structures  with 
low  velocity  zones.  The  phase  with  velocity  about  A  km/s  seen  at 
all  distances  is  a  numerical  artifact. 

a)  Model  1.  An  extreme  inversion  (down  to  A  km/s)  at  2  km  depth. 
Prograde  branch  is  too  strong  at  distances  greater  than  10  km. 


b)  Model  2.  Inversion  at  1.5  km.  Prograde  nergy  is  now  too 
small  at  10  km.  This  depth  was  inconsistent  with  linear 
programming  estimate  of  depth  to  reflector  (>  1.6  km). 

c)  Model  3.  Inversion  at  1.75  km.  Prograde  energy  closer  to 
observations,  caustic  due  to  layer  2-3  boundary  is  better  placed. 

Figure  9.  The  four  seismograms  recorded  nearest  to  the  shot  showing  probable 
reflections . 

Figure  10.  Permeability  log  from  the  IRDP  hole.  Note  the  presence  of  highly 
Impermeable  layers  between  1.2  and  1.7  km  depth.  (From  Johnson 
and  Rummel,  this  volume.) 
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